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Abstract 

An approximate hyperspherical many-body theory has been used to calculate 
the heat capacity and the condensate fraction of a BEC with effective repulsive 
interaction. The effect of interactions has been analysed and compared with the 
non-interacting case. It has been found that the repulsive interaction lowers the 
critical temperature from the value found in the non-interacting case. The dif- 
ference between the critical temperatures increases with the increase in the total 
number of atoms in the trap. 

The number density (ud) of a typical Bose-Einstein condensate (BEC) is re- 
stricted to 10^^ to 10^^ atoms /cm^ so that the average interparticle distance is much 
larger than the range of interatomic interactions. The dilute gas undergoes BEC below a 
certain critical temperature (Tc ~ nano Kelvin), when most of the atoms go to the single 
particle ground state. In this state the momenta of the particles are extremely small and 
the thermal de Broglie wavelength of all particles overlap. The system therefore behaves 
as a single quantum object. At higher temperatures, atoms get distributed into various 
low-lying energy levels. 

It is generally stated that a phase transition occurs during the formation of BEC. 
In contrast with typical classical phase transitions, the origin of this one is the quantum 
mechanical effects of Bose-Einstein statistics. The discontinuity of the heat capacity or 
its temperature-derivative is one of the major manifestations in a phase transition. The 
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heat capacity (Ca) of a system of large but finite number of non-interacting bosons [A) 
in a three dimensional harmonic trap shows a discontinuity at the transition temperature 
(r°), in the semi-classical approximation in which the sum over states is replaced by an in- 
tegral over energy [Il[2]. In this approximation, the chemical potential at temperature T, 
/i(T), is constant (=its maximum value) for T < T° and decreases suddenly for T > T°. 
This approximation is valid for a very large number of particles {A — > oo). Very rapid 
change of Ca{T) near the critical temperature results for the trapped non- interacting 
bose gas, having a large but finite number of particles when sums are numerically evalu- 
ated [3J. The rapid change appears to approach a discontinuity as A ^ oo. Thus for a 
finite system, strictly speaking there is no phase transition, although a rapid change in 
phase occurs near T°, which is the transition temperature for non-interacting bosons in 
the thermodynamic limit. Clearly there is no strict transition temperature for the finite 
system (interacting or non-interacting), a critical temperature (Tc) may be defined as the 
temperature at which Ca attains its maximum In this work we critically examine 
the nature of these quantities for a finite condensate of trapped interacting bosons. 
Interatomic interactions are known to have appreciable effects on the static properties 
of the condensate [IJ. Thus it is important to study the effects of two-body interactions 
on the heat capacity and condensate fraction of the BEC. The most common proce- 
dure is to solve the Gross-Pitaevskii (GP) equation, which is obtained from the mean 
field approach, together with the assumption of a contact two-body interaction, whose 
strength is given by the s-wave scattering length (a^) [T]. A contact interaction is a good 
approximation only in the low density limit. Use of a contact interaction in the mean 
field theory reduces the mean field equation to the GP equation, which is a single second 
order differential equation, non-linear in the condensate wave function. For a finite range 
interaction, the ideal procedure would be to solve the many body Schrodinger equation ab 
initio. Alternatively, one can use an approximate approach like the fully self-consistent 
mean field theory. An exact solution of the Schrodinger equation is impractical for a 
large number (~ 10^) of atoms in the condensate. The essentially exact diffusion Monte 
Carlo (DMC) method [4J is a powerful tool for the many-body problem, but it is rather 
slow and faces difficulties especially for highly excited states of a condensate containing 
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a large number of particles. This is a serious difficulty, since one needs to calculate a 
large number of excited energy levels of the system to obtain thermodynamic quantities. 
In this communication, we have adopted an approximate but ab initio solution of the 
many-body Schrodinger equation, expanding each Faddeev component of the many-body 
wave function in a subset [called potential harmonics (PH)] [5j of the full hyperspherical 
harmonic (HH) basis [6] . The approximation involves disregard of higher-than- two-body 
correlations in the Faddeev component [5], which is well justified in a fairly dilute BEC [7j. 

We have calculated a large number of energy levels (Eni) of a condensate of ^^Rb 
atoms trapped in a spherically symmetric harmonic oscillator potential and using these 
the heat capacity and condensate fraction of the system in the condensed as well as the 
normal phase. Here Eni is the energy in oscillator units of the n^^ radial excitation of the 
Z*'^ surface mode. The energy eigenstates of the system have been calculated using the Po- 
tential Harmonic Expansion Method (PHEM) for trapped interacting bosons ^i8j. This 
technique has been shown to reproduce known results for the static properties [3, El E] 
as also the collapse of attractive condensates [10]. The PHEM was further used to inves- 
tigate the effect of shape dependence of the two-body potential [11] , as also for studying 
the effect of anharmonic traps [T^. These applications have proved that the underlying 
methodology of the PHEM produces reasonable results for the T = properties of the 
BEC. In the present work, the method is extended to investigate thermodynamic quan- 
tities. Convergence of the partition function for T > requires the calculation of a large 
number of energy levels (typically n ~ 400 and / ~ 200). This is very time consuming for 
the essentially exact DMC method. Even for the mean field theory and the CP equation, 
this is a formidable task. By contrast the PHEM is a fairly fast procedure and such a 
calculation is within the realm of feasibility. 

Here, we consider a system of A = -|- 1 identical bosons, each of mass m and 
confined in a trap which is approximated by a spherically symmetric harmonic oscillator 
potential with frequency u. For the static properties, it is assumed that the atomic cloud 
is at zero temperature. The time independent Schrodinger equation is given by, 

„ 4 O „ ,1 J- 



^(f) = E'^{x) (1) 
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where x = {xi, X2, , xa} represents the position coordinates of A particles, V{xi — £,) 

is the pairwise local central two-body interaction between the i^^ and j^^ particles and E' 
is the total energy of the system. The centre of mass motion can be decoupled and the 
Schrodinger equation for relative motion of the system is expressed in terms of Jacobi 
coordinates <^i{i = l, N) > (which are linear combinations of the position coordinates [6]) 



^(ei,6,...,6v) = o (2) 



as 

-|^S-iV| + S-^m^^C^ + S^^.>,V(r-;,) - E 

where E = E' —energy of centre of mass motion and r = [S^^,^?]^/^ is called the hyper- 
radius. The evolution of the system can be studied by following the motion of one point 
in the 3A^ dimensional hyperspace. The polar coordinates of this point are given by a set 
(fijv) of (3A^ — 1) angles. We choose fij = ^^v- For the remaining (A^ — 1) Jacobi vectors 
a hyperradius pij is defined in 3(A^ — 1) dimensional space by pij = [T.^^^^IY^'^ . In the 
PHEM, only two-body correlations are incorporated in the wave function. Higher body 
correlations can be neglected since the gas is very dilute and the probability of a three 
body collision is minimal. So, the wave function \E'(^) can be decomposed into Faddeev 
components 

^iO = ^tJ>^Mr^,,r). (3) 

The Faddeev component ipij describes the (partial) motion of the system when the ij-pair 
interacts, while the remaining {A — 2) particles are simply spectators. The Schrodinger 
equation for the Faddeev component can be written as 
n"^ 1 



2^-»=, ■ . -^=X2 -^ - i^ijinj^r) = -V{rij)J:^^i^f^ijjki{fki,r). (4) 

Summing eq. (4) over all pairs and using eq. (3), we get back the full Schrodinger equation. 
Assumption of two-body correlations alone makes ipij a function of fij and r only [5J and 
hence ipij can be expanded in the complete set of potential harmonics (which are the 
subset of full hyperspherical harmonics, needed for the expansion of V{rij) p]) as 

r) = r^^Sx/t^+i(a,)M^(r), (5) 

where K is the grand orbital quantum number and Qij is the set of all hyperangles 
for the particular choice = Vij. Substitution of eq. (5) in eq. (4) and subsequent 



projection on the PH basis leads to a set of coupled differential equations CDE [71 [8] 
which are then solved using the hyperspherical adiabatic approximation (HAA) [T3l [T^. 
The latter assumes that the hyperradial motion is slow compared to the hyperangular 
motion. This approximation has been shown to be very reliable in atomic and molecular 
cases [13j. The adiabatically separated hyperangular eigenvalue equation is solved (by 
diagonalizing the corresponding potential matrix) to obtain the lowest eigenpotential 
ujQ{r) as a parametric function of r. This is the effective potential for the condensate 
to move as a single quantum entity in the hyperradial space. In the HAA approach, an 
approximate solution of the CDE is obtained by solving a single uncoupled differential 
equation, 

Co(r) = (6) 



m dr"^ dr 
where Co{r) is the condensate wave function in the hyperradial space and {xk,o} is the 
eigen column vector, corresponding to the lowest eigenvalue uJo{r), of the potential matrix 
for a fixed value of r. Ground state in the effective potential well uJo{r) gives the ground 
state energy (Eqq) of the condensate. For the calculation of thermodynamic properties 
using the grand canonical partition function, we need a large number of excitation levels 
in this effective potential well which depends on the orbital angular momentum (/) of 
the system. However for / > 0, computation of the potential matrix element is very time 
consuming. On the other hand, the hyper-centrifugal repulsion term appearing in the 
matrix to be diagonalized, is very large for large A compared to the contribution coming 
from V{fij) for I > 0. Thus the hyper-centrifugal repulsion term contributes most to the 
full matrix [15]. Hence contributions to the off-diagonal matrix elements arising from 
/ > are disregarded for the calculation of Eni- Contributions coming from all terms for 
/ = are properly taken [15]. Finally, the hyperradial equation is solved in the extreme 
adiabatic approximation [13] to calculate the ground and excited energy levels of the 
condensate. 

We perform the calculations for a condensate of ^^Rb atoms with = 2.09 x 10~^ 
o.u.(6.39 X 10"^'' m), which is within the range of values of used in the JILA experi- 
ment [16]. We select only one typical value of the repulsive s— wave scattering length to 



5 



demonstrate our results, since we need to calculate a large number of energy levels, which 
is quite time consuming even by the PHEM. More detailed calculations, particularly those 
for attractive (negative as) condensates will be undertaken later. Although an axially 
symmetric trap (with radial and axial frequencies Ur and Ua respectively) was used in the 
JILA experiment, we assume a spherically symmetric trap of frequency to = (cj^cu^)!, for 
simplicity and to keep our calculations manageable. The interatomic potential is chosen 
to be a realistic one, viz., the van der Waals potential, with a hard core of radius Tc 

V(rij) =00 ,rij<rc 

= ,ri,>r, (7) 

ij 

The value of Cq is known for rubidium atoms [2]. Oscillator units (o.u.) are used in 
our calculations: hu for energy and for length. Value of Cg is 6.489755 x 10"" 

o.u. (4.466 X 10^^^ J.m^). Since the binary collisions occur at extremely low energy, the 
effective atom- atom interaction is specified by the s-wave scattering length a^, which in 
turn depends strongly on Tc. The zero energy two-body Schrodinger equation is solved to 
obtain analytically |2]. The value of Vc is adjusted such that has the experimental 
value. Corresponding two-body wave function is used as a short range correlation function 
for the PH expansion, to enhance its convergence rate [10]. The expansion basis is then 
truncated subject to the condition of convergence of the T = static properties of the 
condensate. Next, a large number of energy levels of the condensate are calculated for 
each of the orbital angular momenta from / = to 200. Calculation of a large number 
of energy levels is very time consuming. Hence for each value of /, a smaller number of 
low-lying levels were calculated directly solving the hyperradial equation. These were 
then least square fitted to a suitable power series expansion. Convergence of such an 
expansion upto the desired accuracy was ascertained. Using this, high-lying levels are 
then obtained by extrapolation. 

The Bose distribution function, f{Eni), is given by 

f (^«') = e/3(Sn.-/.)_i' 
where d = - — — , being the Boltzmann's constant, T is the absolute temperature 



and /i = (T) is the chemical potential. Since the number of bosons {A) is fixed, /i is 
obtained from the constraint [3] 

A = ^Zo^n=odif{Eni), (9) 

where = 2/ + 1 is the degeneracy factor of the l*^ surface mode. The total energy 
E{A,T) of the system is given by 

E{A, T) = ^r=,^'^=Af{Eni)Eni. (10) 

Sums in eqs. (9) and (10) are truncated after achieving convergence upto desired accuracy. 
The heat capacity of the system, C^(T), for fixed particle number {A) is given by 

CAT)^'-^ (11) 
where is given by differentiating eq. (9) with respect to T 



g/z _ T.'^^^T.'^^MEmk - /i)e^(^-^-^)[/(ij;^fc)] 
dT TS~oS-oC?^e/^(^"'-^)[/(^nz)]2 



(13) 



We look for convergence of the chemical potential, as (n, /) sums are truncated in the 
double sum in eq. (9). This v 
eqs. (13) and (12) respectively. 



djji 

double sum in eq. (9). This value of /i(r) is used to calculate and Ca(T), using 



In Fig.l we plot reduced chemical potential, yu/yUo (where yUo is the chemical 
potential at T = 0) as a function of the reduced temperature T/T° (T° is the reference 
critical temperature, ksT^ = 0. 94:1110 A^^^, according to eq. (2.20) of Ref. [2]) for A = 
1000, 2000, 3000,4000 and 5000. Note that in the text book treatment p], /i is taken to 
be equal to the ground state energy of the system for T < T° and it suddenly starts to 
differ for T > T°. In our treatment, since A is relatively small, we evaluate the sums 
over n and I explicitly and ;u(T) is determined from the condition (9) for all T > 0. As a 
consequence /i(T) is a continuous function of T, although /j, remains practically constant 
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Figure 1: (Colour online) Chemical potential calculated by PHEM as a function of tem- 
perature for indicated number of interacting bosons. 



over a wide range of T/T°, upto T/T° ~ 0.8. As T/T° approaches 1, yu/yUo decreases 
rapidly. Also with increasing A the deviation of fi/fio from 1 becomes more sudden. It 
appears that has a sudden change as A ^ oo. 

In Fig. 2, we plot — — as a function of T/T^ for the same number of particles. 



Ak 



B 



The general pattern is similar to the non-interacting case [3]. One notices a sharp change 
just above a critical value (Tc) of T. There is a distinct peak in Ca{T). In 



m 
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B 



the absence of a discontinuity in Ca or its temperature derivative, we follow Ref. [3] to 

OCa 

define the critical temperature (Tc) to be the temperature at which ^„ = 0. It can 



Ca 

be seen from Fig. 2 that T /T^ for which 



Ak 



dT 

is maximum decreases with A, although 



B 



Tc increases with A. Calculated values of Tc are comparable with the experimental 
data [16j. Although a measurement of Tc for such a small number of atoms has not been 
reported, the temperature at which BEC formation was initiated for a larger {A > 10000) 
number of particles in the trap was reported to be about 15nK [16]. The values of 
Tc are listed in table 1, together with T° and T^ (critical temperature for a cloud of 
A non-interacting bosons in an isotropic harmonic trap). It is seen that the effect of 
interaction lowers the critical temperature. This is similar to the result obtained from 
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Figure 2: (Colour online) Heat capacity calculated by PHEM as a function of T/T° 
for indicated number of interacting bosons. A ^ oo indicates infinite number of non- 
interacting bosons. 

the GP equation [1], although the amount of decrease is different. The observation that 
the critical temperature of the interacting gas decreases compared to the non-interacting 
atoms is in conformity with other theoretical and experimental findings p!71 [18] . The 
effective repulsive interaction increases the energy of the system; the system therefore 
has to be cooled to even lower temperatures for all particles to be in the ground state. 
As T increases above T^, most of the atoms get distributed in higher energy levels, with 
a microscopic fraction of atoms left in the ground state. At T = T^, the number of atoms 
left in the ground level is still appreciable for small A - it is denoted by Aq{Tc) and 
presented in the last column in Table 1. Although Aq{Tc) increases with A, the relative 

fraction ^ dec.ea=e= with A 

A 

The variation of critical temperature with the number of bosons in the condensate 
has been presented in Fig. 3. Dependence of T°, and Tc on the number of bosons (A) 
are depicted by curves labelled as 1, 2 and 3 respectively. The interparticle interaction 
is switched off while calculating T^. The system then effectively reduces to A identical, 
non-interacting bosons in an isotropic, harmonic potential, which is identical with the 
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calculation of Ref . [3] . However, both the effects of finite particle number and interparticle 
interactions are included in the calculation of Tc. One notices that {T^ — Tc) increases 
with A. This is intuitively expected since the number of two-body interaction bonds 
increase as A{A — l)/2. 
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Figure 3: (Colour online) Critical temperature as a function of number of Bosons con- 
fined in a spherically symmetric harmonic oscillator trap : 1 - non-interacting Bosons 
in the thermodynamic limit (T^), 2 - finite number of non-interacting Bosons (T/), 
3 - interacting Bosons (Tc) by PHEM. 



Finally, we calculate the condensate fraction (Fc) as 

A-E' ,dif(E„i) 
FciT) = nj^l^^ (14) 

where S' indicates sum over all (n, Z) except the ground state (0,0). We plot it as a 
function of T in Fig. 4 for condensates with 1000, 2000, 3000, 4000, and 5000 particles. 
These plots are again similar to those for the ideal non-interacting case [H [2]. Some 
fluctuations, especially for larger number of particles at lower temperatures, are due 
to numerical errors. In the same figure, we also plot the condensate fraction of non- 
interacting bosons in the thermodynamic limit (indicated hj A ^ oo). While Fc{T) for 
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Figure 4: (Colour online) Condensate fraction calculated by PHEM as a function of 
Tc/T^ for various indicated values of number of interacting bosons. A ^ oo indicates 
infinite number of non-interacting bosons. 

the non-interacting system in the thermodynamic limit reaches zero sharply at T = T^, 

that for the finite interacting system decreases fairly gradually for T > Tc, after a sharp 

drop at ^ Tc- Also as net interaction increases due to increase in A, the curves are shifted 

further to the left, in agreement with the shift of Ca{T) (Fig. 2). In Fig.5, we compare 

the condensate fraction by the PHEM with that obtained using other approaches [Tl [19]. 

The curves are plotted for a particular value of the dimensionless interaction parameter 
fi{N,T = 0) 



V 



-. The parameter 77 is a measure of the ratio of interaction energy and 



the thermal energy pj. The larger the value of 77, the larger is the interaction energy. 
The PHEM result for A = 5000 corresponds to rj = 0.24 and is plotted against the 
reduced temperature T/Tj^, labelled 'PHEM'. For comparison with the mean field local 
density result, we use eq. (122) of Ref. [I] for 77 = 0.24 to plot the curve labelled 'GP'. 
The effect of two-body interaction reduces the condensate fraction appreciably. The con- 
densate fraction in a BEC with interacting bosons has also been calculated using the 
canonical ensemble [19] (eq. (45) of Ref. [19]) and labelled 'canonical' in Fig. 5. In the 
same figure we also plot Fc{T) for non-interacting bosons in the thermodynamic limit 
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which is labelled 'non-int'. The thermodynamical properties calculated using the canon- 
ical ensemble coincides with those obtained using the grand canonical ensemble in the 
thermodynamic limit. However, the effect of finite number of particles has been incorpo- 
rated in the curve marked 'canonical' and therefore differs from the curve obtained using 
the mean field, local density approach in the grand canonical ensemble. Contact inter- 
action has been used for both these approaches. In our method we use the realistic van 
der Waals potential, incorporate two-body correlations in the condensate wave-function 
and compute thermodynamical properties for finite number of interacting bosons. The 
difference of our results from those obtained using the other approaches can be attributed 
to these causes. 




Figure 5: (Colour online) Calculated condensate fraction as a function of T /T^ for the 
dimensionless parameter r] = 0.24, by different indicated methods. 

To summarise, we have calculated the chemical potential, condensate fraction 
and heat capacity of a condensate containing a fixed number of atoms as a function of tem- 
perature (T), using static energy levels calculated by the potential harmonic expansion 
method. A realistic interatomic interaction viz., van der Waals potential (whose short 
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Table 1: Critical temperatures in nK for different values of A and the remaining number 
of bosons in the lowest energy state at T = Tg. 



A 


T 


C 


C 




1000 


4.89 


5.76 


5.18 


102 


2000 


6.04 


7.26 


6.68 


168 


3000 


6.75 


8.31 


7.73 


225 


4000 


7.26 


9.15 


8.51 


268 


5000 


7.65 


9.86 


9.27 


317 



range behaviour is adjusted to give the correct experimental s-wave scattering length) 
is used as the two-body interaction. In this hyperspherical many-body method all the 
two-body correlations are appropriately taken care of, but higher-than-two-body correla- 
tions are disregarded, which is justified for the dilute condensate. Calculated Ca shows 
a gradual increase with T, until it reaches a maximum and falls rapidly near the critical 
temperature (Tc). There is no discontinuity either in Ca or its temperature derivative 
as functions of T. The sharpness of the sudden fall increases with A. This is similar to 
the non-interacting inhomogeneous BEC where Ca appears to have a discontinuity at 
T — as A ^ oo (7° is the critical temperature in the thermodynamic limit) . Beyond 
Tc, Ca/A approaches the ideal value Sks- We notice that Tc increases gradually with 
A, which is also seen for non-interacting atoms. We find that critical temperature for 
interacing atoms is lower than that of non-interacting atoms, which agrees with intuitive 
expectations. 
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